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SUPERFLUID TURBULENCE AND PULSAR GLITCH STATISTICS 
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ABSTRACT 

Experimental evidence is reviewed for the existence of superfiuid turbulence in a differentially rotating, 
spherical shell at high Reynolds numbers (Re }t 10 3 ), such as the outer core of a neutron star. It is 
shown that torque variability increases with Re, suggesting that glitch activity in radio pulsars may 
be a function of Re as well. The Re distribution of the 67 glitching radio pulsars with characteristic 
ages r c < 10 6 yr is constructed from radio timing data and cooling curves and compared with the 
Re distribution of all 348 known pulsars with r c < 10 6 yr. The two distributions are different, with 
a Kolmogorov-Smirnov probability > 1 — 3.9 x 1CU 3 . The conclusion holds for (modified) Urea and 
nonstandard cooling, and for Newtonian and superfiuid viscosities. 

Subject headings: dense matter — hydrodynamics — stars: interior — stars: neutron — stars: rotation 



1. INTRODUCTION 

The known number of glitching radio pulsars has 
quadrupled in the last four years, as has the to- 
tal number of glitches, in the wake of the Parkes 
Multibcam Su rvey and imp roved multifrequency tim- 
ing solutions (jPeraltal 120061 ). Out of a total popu- 
lation of 1730 pulsars, 101 have experienced a total 
of 287 glitches (IWang et all l2000t iKrawczvk et all l2003t 



iHobbs et all l2004t Uanssen fc Stappersl I2006D . The en 



present numerical simulations of a 1 So-paired neutron su- 
perfiuid in the outer core of a neutron star as the crust 
spins down electromagnetically. We show that the torque 
is steady for Re 10 3 and unsteady (e.g., oscillatory) oth- 
erwise. In we combine standard cooling curves with 
measurements of the characteristic ages of radio pulsars 
to compute the temperature, viscosity, and hence Re of 
the star. We show that the Re distributions of those pul- 
sars that glitch, and those that do not, are different — a 
rare instance of regularity in the glitch phenomenon. The 
implications for the existence of superfiuid turbulence in 
neutron stars are weighed critically in fJU 

2. HYDRODYNAMIC SUPERFLUID TURBULENCE 

Large-scale, time-dependent, meridional circulation, 
driven by Ekman pumping, is a generic feature of viscous 
and superfiuid flow in spherical Couette geometry (i.e., a 
differe ntially rotating, spherical shell); see Uunk fc Egbersl 
(|2000f) for a review. The complexity of the global circula- 
tion pattern increases with Re. When Re is large enough, 
the f low becomes nonax i symm etric and ultimately chaotic, 
e.g., iNakabavashi et all ([20021) . 

Figure [1] compares spherical Couette flow of a two- 
component, Hall-Vincn-Bekarevich-Khalatnikov (HVBK) 
superfiuid at low and high Re. Meridional streamlines of 
the normal (i.e., torque generating) component are plot- 
ted for Re = 3 x 10 2 and 3 x 10 4 in Figures la and lb 
respectively. The flow is computed by integrating the 
HVBK equations (jHills fc Roberts! [l977l ) n umerically us- 
ing a pseudospectral collocation method (jCanuto et al.l 
|1988|) with spectral filtering to su ppress high spatial 
frequen c ies, e specially at high Re (|Peralta et alj 120051 : 
iPeraltal [2006). We adopt no-slip and no-penetration 
boundary conditions for the normal component (veloc- 
ity v„), perfect slip for the superfiuid component (i.e., 
no pinning, for numerical stability; velocity v s ), and a 
Stokes flow (with v„ = v s ) initially. The mutual friction 
force is of the isotropic . Gorter-Mellink form ( oc vf ls v ns , 
with v„ s = v„ — v s ) (|Gorter fc Mellinkl Il949f ). because 
the meridional counterflow in Figures [T^ and [T}) excites 
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larged data set invites an updated statistical study. Pre- 
vious analyses identified certain empirical trends, e.g., 
glitch acti yity and healing parameter versus character- 
istic age (IShemar fc Lvnel ll996t lUrama fc Okekd 119991: 
I Wane: et alJl200Cfl. and glitch amplitude versus recurrence 
time ( Middleditch ct al. 2006), while refuting certain theo- 
retically predicted correlations, e.g., glitch amplitude ver- 
sus spin period P, period derivative P, an d total spin down 
since the last glitch ("r eservoir model") |Sh emar fc Lvnd 
[19961 IWang etalll2000h . 

On the theoretical front, the global pattern of superfiuid 
circulation in the outer core of a neutron sta r was recently 
computed n umerically for the first time (jPeralta et alj 
I2005L l2006al fbl). In general, a fluid contained in a differ- 
entially rotating, spherical shell flows unsteadily under a 
wide range of conditions, especiall y in thick shells and a t 
high Reynolds numbers Re > 10 3 (|Junk fc Egbersl l2000f) . 
Such unsteady flow is observed in numerical and labora- 
tory experiments with visco us fluids, e.g., transitions be - 
tween Taylor vortex states dMarcus fc Tuckermanl Il987| ) , 
Taylor vortex oscillations (iHollerbacbl 1998h . relaminar- 
ization ([Nakabavashi et al.l I2005T T and nonaxisymmetric 
modes on the route to h igh-i?e turbulence, like Taylor- 
Gortler vortices (jLj|[2004). Time-dependent flow is also 
observed in superfluids like helium II, e.g., torque "steps" 
and quasiperiodi c torque oscillations in Coue tte and spin- 
up experiments (jTsakadze fc Tsakadzelll980h . 

In this Letter, we combine theory and observation to 
test whether superfiuid tur bulence exists in n eutron stars, 
as originally proposed by iGreensteTnl (|l97Gj ). In SJ2 we 
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Kelvin waves via the Donnelly-Glaberson instability, con- 
verting the rectilinea r vortex array into a vortex tangle 
(ISwanson et alj Il983t iTsubota et all 120031: IPeralta et all 
2006b). The superfluid component also feels a vortex ten- 
sion force oc ijj s x (V x u> s ), with oj s = V x v s . 

The meridional circulation pattern in Figure [Tt> is more 
complicated than in Figure [Tb; seven cells are visible for 
Re = 3 x 10 4 , compared to only one cell for Re = 3 x 10 2 . 
The azimuthally aver aged differential torque dN z /d(cos9) 
(|Peralta e t al. 2006b) is a monotonic function of latitude 
for Re — 3 x 10 2 but peaks near the (moving) secondary 
vortices (at 9 w 30°) in Figure [Tb for Re = 3 x 10 4 . 

The torque is variable at high Re, exhibiting persistent 
oscillations under steady differential rotation , whose am- 
plitude depends on Re (|Peralta et alj [20051 ). Figure [Tb 
displays the torque on the outer sphere (stellar crust) as 
a function of time. For t > 20 fi -1 , with £1 = 2tt/P, 
after initial transients die away, we see that the torque 
is asymptotically steady for Re = 3 x 10 2 (solid curve) 
and oscillatory for Re = 3 x 10 4 (dashed curve). The 
variable torque leads naturally to rotational irregularities. 
The variability time-scale (~ Q^ 1 ) is too fast to account 
for pulsar timing noise directly. However, the torque os- 
cillates quasiperiodically, not periodically, so Q wanders 
stochastically when integrated over long times [c.f. micro- 
jumps and osc i llations superimpose d on a random walk; 
iBovnton et ail (| 19721) ; iHobbi (12002ft ]. 

Of course, we wish to investigate the flow at 
Re ^ 10 9 , the realistic R e regime for radio pulsars 
(Ma strano fc Melatosl 120051 ). Unfortunately, our simula- 
tions cannot handle this. Even for Re — 10 5 , and with 
heavy spectral filtering, the flow is unresolved numeri- 
cally; the spectral coefficients do not dec line monoton i- 
cally with polynomial order as they should ( Peralta 2006) . 
From laboratory experiments, however, it is clear what 
to expect qualitatively as Re increases: the flow passes 
through a sequence of bifurcations to nonaxisymmetric 
state s like herringbone wav es and Stuart vortices at Re ~ 
10 6 ijJunk fc Egbersl I2000D . before developing fully into 
scale-free, Kolmogorov-like turbulence at R e ^ 10 7 , which 
itself is not always isotropic at small scales (jDouadv et alj 
1991). Turbulent flow in He II at Re ^ 10 5 often contains 
a vortex tangle (fBarcngh i et al.iri997ft . Moreover, in ex- 
periments with He II, eddies in the normal and superfluid 
components match at scales large r than the average vortex 
separation (|Barenghi et al.l [19971 ). implying that, macro- 
scopically, high-i?e superfluid turbulence behaves similarly 
to high-i?e classical turbulence. Hence, in neutron stars 
(Re ^ 10 9 ), the vortex tangle in the outer core may match 
its vorticity to the normal component. The vortex tension 
modifies this picture somewhat by contributing "stiffness" 
in regions where u> s builds up (e.g., at the smallest scales 
in the Kolmogorov cascade). The largest eddies in fully 
developed turbulence, where the kinetic energy mainly re- 
sides, are nonaxisymmetric and move "jerkily" , so the net 
torque fluctuates substantially. 

3. Re DISTRIBUTION OF GLITCHING PULSARS 

Given that the Reynolds number is a key factor govern- 
ing the variability of the global superfluid flow in a neutron 
star, and hence the star's rotation, it is interesting to test 
whether the amplitude and rate of incidence of rotational 



irregularities depends on it. From the definition 

Re = R 2 VL/v n , (1) 
where R is the stellar radius and v n is the kinematic 
viscosity of the normal component, it is clear that Re 
can be measured in principle, if we know how v n de- 
pends on the density p and temperature T in the outer 
core. Using the Newtonia n viscosity formula derived by 
ICutler fc Lindblom (1987), due to neutron- neutron scat- 
tering, we find 

Re = 1.8 x 10 11 (T/10 8 K) 2 (f7/10 2 rads- 1 ), (2) 
with p = 2.8 x 10 12 g cm~ 3 (normal component) and R = 
10 6 cm. For electron-electron scatterin g in a superfluid, 
the m ultiplicative factor is 6.0 x 10 10 (jAndersson et al.l 
2005). In writing equations flTJ) and ([2]), we imagine a 
model where the outer core extends over the density range 
(2-4) xlO 14 g cm~ 3 , where the superfluid behaves hy- 
drodynamically (i.e., weak bulk pinning) and isotropically 
( x So-paired). We assume that the normal component con- 
stitutes 1 % of the total density, to account for conden- 
sate depletion by nonideal effects. We do not treat strat- 
ification in the simulations to keep them tractable, al- 
though stratification is known t o be important in sup- 
pressing meridional circulation (|Abnev fc Epstein] [l996; 
IPeralta et ai1l2006bft . 

The core temperature T is related to the surface tem- 
pe rature T s , e.g., via t he tw o-zone, heat-blanket model 
of iGudmundsson et all (|1982h . which gives (T/10 8 K) = 
1.29(T S /10 6 K) 8 . Surface temperatures have been mea- 
sured for 11 neutron stars, by fitting their thermal X- 
ray emission to hydrogen or heayy-elem ent atmospheres 
(|Larson fc Linkl l2002t iPage et alTi2004ft . However, the 
sample is too small to analyse statistically; only three 
of these objects exhibit glitches, for example. Instead, 
we are forced to estimate T s from the characteristic age 
t c = P/(2P), combined with theoretica l cooling c urves, 
calibrated against the 11 objects above (|Pagdfl9 9ll). We 
consider three cooling mechanisms: (I) standard neu- 
trino cooling, via th e Urea and modified Urea processes 
(Hof fberg et al.|[i970T ). (II) nonsta ndard (fast) cooling due 
to pa ired neutron superfluidity (|Amundsen fc Ostgaardl 
1985j), and (III) nonstandard cooling with the superfluid 
energy g ap fitted empirically to observations of therma l 
emission (IGnedin et aLlll99l : iLevenfish fc Yakovlevlli996l) . 

Figure [2b displays the raw Re distributions of the 67 
glitching pulsars with r c < 10 6 yr (green histogram) and 
all 348 known pulsars with r c < 10 6 yr (red histogram), 
including the 67 glitchers. The objects are sorted into ten 
bins of equal logarithmic width 0.37 dex in the interval 
9.4 < log 10 (i?e) < 13.1. Standard cooling (I) and Newto- 
nian viscosity formulas are applied to both samples. We 
restrict attention to object s with r e < 10 6 yr, because the 
cooling curves compiled by I Page! ([1998) drop away precip- 
itiously above this age, as neutrino cooling gives way to 
photon cooling (which depends sensitively on the uncer- 
tain composition of the stellar atmosphere). 

The histogram in Figure [2b is redrawn for nonstandard 
cooling (II and III) in Figures [5h> and [2b, with bin widths 
of 0.30 dex and 0.96 dex, respectively. Identical (after 
rescaling Re) distributions are obtained if the superfluid 
electron-electron scattering viscosity is used instead. 

The Re distribution of the glitching pulsars differs 
clearly from that of the total population. One sees im- 
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mediately, from Figure [5^,, that most of the glitching pul- 
sars have 10 11 ^ Re ^ 10 12 , whereas the pulsar population 
peaks at 10 10 £ Re <, 10 . To quantify this, we perform a 
Kolmogorov-Smirnov (K-S) test on the cumulative distri- 
butions constructed from Figure^ (to circumvent binning 
bias). The K-S statistic is D = 0.3, yielding a probability 
p = 1 — 1.6 x 10 -5 that the two populations are drawn 
from different distributions. Similar conclusions follow for 
nonstandard cooling of types II and III, and for superfluid 
electron-electron scattering viscosity (p does not change 
when Re is rescaled multiplicatively). The associated IC- 
S-probabilities are quoted in Table [1] 

Is Figure [2] just a restatement of the well known trend 
that glitching pulsars tend to be young? Partly, but not 
entirely. The period distribution of glitching pulsars dif- 
fers clearly from that of the total population, with a K- 
S probability p = 0.994. The maximum separation be- 
tween the cumulative P distributions occurs at an inter- 
mediate period P — 0.32 s, ruling out a pure age effect. 
[Observational studies of kick velocities and pulsars in 
supernova remnants suggest that the distribution of the 
birth periods is quite flat in the range 0.05 — .5 s; se e 
iFaucher-Giguere fe Kaspl (|2006l ). iNg fc Romanj (|2007f ). 
and references therein.] It appears that P and r c control 
glitch behaviour independently (perhaps, but not neces- 
sarily, through Re). Interestingly, glitching pulsars are 
distributed differently in r c and Re, peaking at t c ~ 10 6 
yr (oldest objects) and Re ~ 10 11 (intermediate objects) 
respectively. 

Given that Re distributions of glitching and nonglitch- 
ing pulsars are different, it is natural to ask whether the 
rate and amplitude of glitch activity are simple functions 
of Re. There are several possible measures of glitch activ- 
ity in an individual pulsar. We pick three: (i) the activity 
parameter A g — (N g ACl g )/(t g fl), where N g is the number 
of glitches observed over a total observation time t g (taken 
to be the time since discovery, in the absence of complete 
data), a nd A£l a is the accumula ted change in due to 
glitches (jMcKenna fc Lvnel [l990h ; (ii) An g /fl; and (hi) 
the mean recurrence time t r between glitches, normalized 
by t c (in 11 out of 67 objects, this is an upper limit, as 
they have glitched once only). Interestingly, none of these 
measures show a clear trend with Re, leading to scatter 
plots. This suggests that Re controls the threshold for 
glitch activity, rather than its rate and amplitude, which 
is qualitatively consistent with ij^l For example, one can 
imagine a Re threshold (when turbulence sets in), which 
must be exceeded for glitches to occur at all, and a sep- 
arate "rate" threshold (perhaps controlled by the Rossby 
number, proportional to the crust-core shear), which de- 
termines how often the pulsar glitches per unit observa- 
tion time (and hence A g ). Likewise, t r /r c and AQ g /£l are 
related physicallly to the rate at which crust-core shear 
builds up, which is independent of Re. Also, A g is not di- 
mensionless, so it is not surprising that factors other than 
Re affect it. Note that Re decreases with r c b ut varies in- 
depen dently with Q, so the A g -r c correlation (|Wang et al.l 
2000) is washed out in an A g -Re plot. 

4. DISCUSSION 

A superfluid (or, indeed, a Newtonian fluid) contained 
in a differentially rotating, spherical shell circulates merid- 



ionally and exerts a time-dependent viscous torque on the 
outer shell. The vigour of the circulation, and the vari- 
ability of the torque, increase with Re; the flow is steady 
for Re ;$ 10 2 and chaotic for Re ^ 10 5 . Spherical Couette 
flow is an idealised model of the superfluid outer core of a 
neutron star, in which the outer sphere is the stellar crust, 
which spins down electromagnetically. It is therefore in- 
teresting to test whether the incidence of pulsar rotational 
irregularities like glitches depends on Re. We show in JJ] 
that this is indeed the case: the distribution of glitching 
and nonglitching pulsars with t c < 10 6 yr is different, with 
K-S probability > 1 — 3.9 x 10~ 3 . This new fact joins the 
A g -versus-T c correlation as one the few observed regulari- 
ties in glitch phenomenology. 

An immediate worry is that the distributions in Fig- 
ure [5] are contaminated by some observational selection 
effect. Our estimates of Re are uncertain in two ways. 
First, t c is the characteristic spin-down age of the pulsar, 
not its true age. For example, five young pulsars with 
t c < 1.7 x 10 3 yr have braking index n < 3, whereas 
the formula r e = P/(2 P) assumes n = 3 (jMelatos! Il997l; 
iLivingstone et al.ll2006D . Second, the theoretical cooling 
curves depend sensitively on the superfluid energy gap in 
the outer core, which is poorly known (although admit- 
tedly calibrated by 11 pulsars with observed thermal emis- 
sion). However, the above uncertainties do not affect the 
conclusion that the Re distributions of the glitching and 
nonglitching populations in Figure O are different, because 
T and hence Re are calculated in the same way for both 
samples. Indeed, it is hard to imagine that our ability to 
detect glitches in radio timing data is a function of Re. 

Physically, the conclusions from Figure [5] are at once 
natural yet surprising. It is well known that torque vari- 
ability increases with Re in a differentially rotating su- 
perfluid in laboratory and numerical experiments ( - 
However, it is strange that pulsars behave differently at 
Re ?i 10 10 (few glitches) and Re ^ 10 11 (many glitches). 
Naively, Kolmogorov turbulence should be scale- free, fully 
developed, and indistinguishable at these Reynolds num- 
bers. There are at least two ways to explain this. First, 
processes involving m i crosc opic superfluid turbulence, e.g., 
pinning (lAlpar et alJ Il984f ) and vortex tangle formation 
(|Peralta et al.l l2006bf ). may depend sensitively on Re in 
the range 10 10 < Re < 10 11 . Second, it is possible 
that theoretical estimates of Re (jAndersson et all 12005: 
iMastrano fc Melatosl [2005h are ~ 10 6 times too high. If 
so, this reduces the peak of the glitching pulsar distri- 
bution to Re 10 5 , exactly where spherical Couette flow 
breaks up into nonaxisymmetric modes before becoming 
chaotic (iJH). We have no reason to doubt the v n carefully 
derived in the literature. However, if the flow is turbu- 
lent, Reynolds stresses may completely overwhelm viscous 
stresses, such that v n is replaced by v n + A, where A satis- 
fies Advi/dxj w (SviSvj), Vi is the mean velocity, and Sv j 
is the fluctuating velocity (zero average ) (lPedloskvlll982T) . 
Typically, we have Ajv n ~ (Sv/v) 2 Re 3> 1. Turbulent 
Reynolds stresses drastically reduce Re in a range of ge o- 
physical and astrophysical applications (|Pedloskvlll982f ). 

Plainly, there exist several high-iie pulsars that do not 
glitch. Observational selection effects aside, there exist 
several theoretical mechanisms that may be responsible. 
First, for spherical Couette flow in the range 10 3 < Re < 
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10 7 , there are sizable intervals of Re where the flow is 
almost la minar, interposed between interv als where it is 
turbulent (jNakabavashi fc Tsuchidal I2005D . Second, for 
Re > 10 7 , where the turbulence is fully developed (Kol- 
mogorov) and "isotropic" when averaged over a turn-over 
time, the torque on the star is steadier than for organ- 
ised modes (e.g., Taylor-Gortler vortices) near the onset 
of turbulence {Re ~ 10 5 ), which are not isotropic on aver- 
age. Third, if atmospheres and planetary interiors are any 
guide, the turbulent Reynolds stresses depend sensitively 
on what nonlinear modes are excited in the turbulence, 
and the effective Re may vary greatly (and unpredictably) 
from object to object, in ways not captured by 1[2"]), Inci- 
dentally, none of the pulsars in our sample have low Re, 
even after adjusting for turbulent Reynolds stresses; the 
minimum is Re ~ 10 3 . Hence there are no examples of 
low-i?e objects that glitch unexpectedly. 

In closing, we emphasize that the results of this paper 
do not prove that superfluid turbulence exists in a neutron 
star, let alone that it controls glitch behaviour. However, 
considerable effort has been expended to identify patterns 



in glitch behaviour, with limited success. The undeniable 
empirical fact that the Re distributions in Figure [2] are 
different is therefore intriguing, especially because Re is 
a dimensionless quantity with deep hydrodynamical signif- 
icance, and because the K-S comparison is robust with 
respect to how T, v n , and hence Re are estimated obser- 
vationally. One hopes it will offer an enlightening clue to 
glitch physics, whether or not superfluid turbulence turns 
out to play a central role. We also hope that the results 
presented here will stimulate the ongoing observational 
campaign to measure pulsar temperatures. 
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- 6.4 x 10~ b 


II 
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-3.9 x 1(T 3 


III 


1 


-2.5 x 10" 3 



Table 1 

k-s probabilities p for newtonian viscosity and three cooling mechanisms. 




Fig. 1. — Snapshot at t = 50f2 — 1 of meridional streamlines (computed by integrating the in-plane velocity vector) for the normal fluid 
component in spherical Couette flow in the outer core of a neutron star for (a) Re = 3 X 10 2 and (b) Re = 3 X 10 4 . (c) Evolution of the outer 
torque on the crust (in units of pi? 5 f! 2 ) for Re = 3 X 10 2 (solid curve) and Re = 3 X 10 4 (dashed curve). The simulation parameters are 
mutually ordered as in the outer core of a neutron star but exaggerated for computational tractability, viz., crust-core angular velocity shear 
AO/fi = 0.3, dimensionless shell thickness 8 = 0.5, stiffness parameter u s = 10~ 5 iJ 2 f! u n , Gorter-Mellink constant A' = 5.8 X 10 — 2 , and 
superfluid density fraction ps/p = 0.99. We follow exactly the notation and definitions in Pcralta ct al. (2005). 




Fig. 2. — Reynolds number distribution for the 67 glitchers (green curve) and all 348 pulsars (red curve, including the 67 glitchers) with 
characteristic age t c < 10 6 yr, using (a) standard neutrino cooling, (b) nonstandard (fast) cooling due to paired neutron superfluidity, and (c) 
nonstandard cooling with the superfluid energy gap fitted empirically to observations of thermal emission. The data are divided into 10 bins 
of equal logarithmic width (a) 0.37, (b) 0.30, and (c) 0.96. The percentage of objects that glitch in each bin, for standard neutrino cooling, 
is recorded above the histogram in (a). 



